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ON THE ORIGIN OF THE GLOBAL SCHMIDT LAW OF STAR FORMATION 
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ABSTRACT 

One of the most puzzling properties of observed galaxies is the universality of the empirical correlation between 
the star formation rate and average gas surface density on kiloparsec scales (the Schmidt law). In this study I 
present results of self-consistent cosmological simulations of high-redshift galaxy formation that reproduce the 
Schmidt law naturally, without assuming it, and provide some clues to this puzzle. The simulations incorporate 
the main physical processes critical to various aspects of galaxy formation and have a dynamic range high enough 
to identify individual star forming regions. The results indicate that the global Schmidt law is a manifestation of 
the overall density distribution of the interstellar medium (ISM). In particular, the density probability distribution 
function (PDF) in the simulated disks is similar to that observed in recent state-of-the-art modeling of the turbulent 
ISM and has a well-defined generic shape. The shape of the PDF in a given region of the disk depends on the 
local average surface density S g . The dependence is such that the fraction of gas mass in the high-density tail of 
the distribution scales as S"" 1 with n « 1.4, which gives rise to the Schmidt-like correlation. The high-density tail 
of the PDF is remarkably insensitive to the inclusion of feedback and details of the cooling and heating processes. 
This indicates that the global star formation rate is determined by the supersonic turbulence driven by gravitational 
instabilities on large scales, rather than stellar feedback or thermal instabilities on small scales. 

Subject headings: cosmology: theory-galaxies: evolution-galaxies: formation-ISM: kinematics and dynamics - 
ISM: structure - stars: formation-methods: numerical 



1. INTRODUCTION 

Formation of galaxies is a complicated process involving a 
variety of physical phenomena on a large range of spatial and 
temporal scales. Although hierarchical structure formation sce- 
nario proved to be a very successful framework for interpreting 
a wide variety of observational data, much remains to be under- 
stood about the processes that shape properties of galaxies. In 
particular, star formation is one of the most important, yet very 
poorly understood processes. 

On small scales (< 10-100 pc), star formation appears to be 
a complicated and stochastic process. Understanding how stars 
and stellar clusters form in individual molecular clouds is still 
sketchy, even in the best studied cases (e.g., Hartmann 2002). It 
is therefore remarkable that on kiloparsec scales star formation 
exhibits surprising regularity. Schmidt (1959) first argued that 
the star formation rate (SFR) in galaxies is a power law function 
of average gas density. This simple 'Schmidt' law 

S SFR «E" g , (1) 
where Ssfr and E g are the surface densities of young stars and 
gas averaged on a large (~ kpc) scale, was found to hold for a 
wide range of galaxy types, star formation rates, and gas den- 
sities (Kennicutt 1983, 1989, 1998a,b). The star formation his- 
tories of the Local Group dwarfs indicate that the Schmidt law 
applied at high redshifts (Gnedin 2000). The power-law index 
of the correlation falls in the range n ~ 1.3 — 1.6, depending on 
the tracers used and the scales considered (Kennicutt 1998a). 
In addition, star formation rate in galaxies declines sharply be- 
low a critical threshold surface density of ~ 5 — 10 M Q pc~ 2 
(e.g., Kennicutt 1989; Martin & Kennicutt 2001, and references 
therein). This density is thought to be associated with large- 
scale gravitational and shearing stability thresholds, although 
details of the mechanism are still uncertain. Nevertheless, at 
high gas densities the form of the Schmidt law appears to be 
remarkably consistent from galaxy to galaxy, both in terms of 
slope and absolute efficiency (i.e., normalization of the rela- 
tion). 

Despite the apparent simplicity, a unique explanation for the 
Schmidt law proved to be elusive. To a large extent, the prob- 



lem is that the Schmidt law can be explained by any process in 
which gas consumption time depends on the local average dy- 
namical time (oc p~ 1 ' 2 ). A wide variety of possible models with 
such scaling have been proposed in the last three decades (see 
Elmegreen 2002, for a review). Nevertheless, the universality 
of the star formation efficiency and the slope of the Schmidt law 
in the face of complexity of star formation on small scales re- 
mains a puzzle. In this Letter I present results of self-consistent 
cosmological simulations of high-redshift galaxy formation that 
reproduce the Schmidt law naturally, without assuming it, and 
provide some clues to this puzzle. 

2. IMPLEMENTATION OF STAR FORMATION 

The simulations presented in this paper were performed us- 
ing the Eulerian gasdynamics+iV-body Adaptive Refinement Tree 
(ART) code. The code is based on the cell-based approach 
to adaptive mesh refinement (AMR) developed by Khokhlov 
(1998). The algorithm uses a combination of multi-level particle- 
mesh (Kravtsov et al. 1997; Kravtsov 1999) and shock-capturing 
Eulerian methods (van Leer 1979; Colella & Glaz 1985) to fol- 
low the evolution of the DM and gas, respectively. High dy- 
namic range is achieved by applying adaptive mesh refinement 
to both gasdynamics and gravity calculations. 

Several physical processes critical to various aspects of galaxy 
formation are implemented in the code: star formation, metal 
enrichment and thermal feedback due to the supernovae type 
II and type la (SNII/Ia), self-consistent advection of metals, 
metallicity- and density-dependent cooling and UV heating due 
to cosmological ionizing background using cooling and heating 
rates tabulated for the temperature range 10 2 < T < 10 9 K and 
a grid of densities, metallicities, and UV intensities using the 
Cloudy code (ver. 96b4, Ferland et al. 1998). The cooling 
and heating rates take into account Compton heating/cooling 
of plasma, UV heating, atomic and molecular cooling. The 
detailed implementation of these processes will be described 
elsewhere (Kravtsov 2003). As I will show below, the results 
presented in this paper are not sensitive to the details of these 
processes. Here I will focus on the implementation of star for- 
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mation, the process which is the subject of this study. 

In numerical simulations star formation is usually assumed 
to occur in certain "star forming" regions identified using some 
specific criteria. The gas in such regions is then converted 
into stars on a characteristic gas consumption time scale, r»: 
p* = Pg/r*. Observationally, the efficiency of star formation on 
small (< 100 pc) scales in the densest regions of the ISM is sig- 
nificantly higher than the low global efficiency indicated by the 
Schmidt law (e.g., Elmegreen 2002). This means that the value 
of t* and, perhaps more importantly, its density dependence 
in simulations should be different depending on whether star 
forming regions are resolved. Clearly, the use of the Schmidt 
law is observationally motivated only on scales > 1 kpc. Al- 
though implementations vary, galaxy formation simulations usu- 
ally adopt gas consumption time that scales as t* oc max(/ coo i,fdynX 
where f coo i and fd yn are the local cooling and dynamical time, 
respectively. In even moderately dense regions, f coo i <C fd yn and 
t* oc fdyn oc /T 1 / 2 , which results in the Schmidt-like star forma- 
tion law: p* oc /9g 5 . 

In this study, a "minimal" star formation prescription was 
used. Namely, the stars were formed with a constant r* so that 
p* oc p g . This "constant efficiency" model on the scale of star 
forming regions is well motivated by observations (e.g., Young 
et al. 1996; Wong & Blitz 2002). The star formation was al- 
lowed to take place only in the coldest and densest regions, 
T < Tsf and p g > psf, but no other criteria (like the collapse 
condition V • v < 0) were imposed. I used r* = 4 Gyrs, 7sf = 
9000 K, and p SF = 1 .64 M Q pc" 3 or atomic hydrogen number 
density of «h = 50 cm -3 . The adopted values of 7sf and psf are 
quite different from the typical temperatures and density of star 
forming molecular cores: T < 30-50 K and «h ^ 10 4 cm" 3 . 
They are, however, more appropriate to identify star forming 
regions on ~ 100 pc scales which are resolved in the presented 
simulations. In practice, T$p is not relevant because most of the 
gas at p > psF is at temperatures of just a few hundred degrees 
Kelvin. The adopted gas consumption time scale t* is derived 
from the observed normalization of the Schmidt law. 

Algorithmically, star formation events are assumed to occur 
once every global time step Afo < 10 7 yrs, the value close to 
the observed time scales (e.g., Hartmann 2002). Collisionless 
stellar particles with mass m* = p* Afo are formed in every un- 
split mesh cell during star formation events. The mass of stellar 
particles is restricted to be larger than 10 4 M but smaller than 
2/3 of the gas mass contained in the parent cell. This is done 
in order to keep the number of stellar particles computationally 
tractable while avoiding a sudden dramatic decrease in the local 
gas density. 

Once formed, each stellar particle is treated as a single-age 
stellar population and its feedback on the surrounding gas is 
implemented accordingly. I use the word feedback in a broad 
sense to include injection of energy and heavy elements (met- 
als) via stellar winds and supernovae and secular mass loss. 
Specifically, in the simulations analyzed here, I assumed that 
the stellar initial mass function (IMF) is described by the Miller 
& Scalo (1979) functional form with stellar masses in the range 
0. 1 - 100 M . All stars with M»>8M deposit 2 x 10 51 ergs 
of thermal energy and a mass fjM* of heavy elements in their 
parent cell (no delay of cooling was introduced in these cells). 
The metal fraction is fz = min(0.2, 0.OlM^-0.06), which crudely 
approximates the results of Woosley & Weaver (1995). In ad- 
dition, the stellar particles return a fraction of their mass and 
metals to the surrounding gas at a secular rate mi oss = m* Q(f- 



fbirth + To) -1 with Co = 0.05 and To = 5 Myr (Jungwiert et al. 
2001). The code also accounts for SNIa feedback assuming 
a rate that slowly increases with time and broadly peaks at the 
population age of 1 Gyr. I assume that a fraction of 5 x 10" 3 of 
mass in stars between 3 and 8 M Q explodes as SNIa over the 
entire population history and each SNIa dumps 2 x 10 51 ergs 
of thermal energy and ejects 1 .3 M Q of metals into parent cell. 
For the assumed IMF, 75 SNII (instantly) and 1 1 SNIa (over 
several billion years) are produced by a 10 4 M Q stellar particle. 

3. NUMERICAL SIMULATIONS 

In this paper I present two simulations of the early (z > 4) 
stages of evolution for a galaxy of typical mass: « 10 12 /i _1 M Q 
at z = 0. At the analyzed epochs, the galaxy has already built 
up a significant portion of its final mass: 1.3 x 10 10 /i _1 M Q at 
Z = 9 and 2 x I0 ll h~ l M Q at z = 4. The evolution is started from 
a random realization of a gaussian density field at z = 50 in a 
periodic box of 6/T 1 Mpc with an appropriate power spectrum 
and is followed assuming flat ACDM model: f^o = 1 - Qa = 
0.3, fl h = 0.043, h = Ho/100 = 0.7, n s = 1, and cr 8 = 0.9. The 
parameters have their usual meaning and are consistent with 
recent cosmological constraints. 

To increase the mass resolution, a low-resolution simulation 
was run first and a galactic-mass halo was selected. A lagrangian 
region corresponding to five virial radii of the object at z = 
was then identified at z = 50 and re-sampled with additional 
small-scale waves (Klypin et al. 2001). The total number of DM 
particles in the high-resolution lagrangian region is 2.64 x 10 6 
and their mass is moM = 9.18 x lOV M n . Outside the high- 
resolution region the matter distribution was sampled with « 
3 x 10 5 higher mass particles. 

As the matter distribution evolves, the code adaptively and 
recursively refines the mesh in high-density regions. In the sim- 
ulations I present, two main refinement criteria were used: gas 
and DM mass in each cell. The code used a uniform 64 3 grid 
to cover the entire computational box. The lagrangian region, 
however, was always unconditionally refined to the third re- 
finement level, corresponding to the effective grid size of 512 3 . 
Beyond the third level, a mesh cell was tagged for refinement if 
its gas or DM mass exceeded 0.125 and 0.0625 times the mean 
mass expected for the average density in each component in the 
zeroth level (i.e., uniform grid) cell, respectively. The refine- 
ment thus follows the collapse of 1.2 x 10 6 /i _1 M Q (gas) and 
3.7 x 10 6 /r' M Q (DM) mass elements in a quasi-lagrangian 
fashion. 

The maximum allowed refinement level f max was set to nine. 
A total of :=y 1.1 x 10 7 mesh cells was used at z = 4 with « 
2.5 x 10 5 of them at refinement levels of 8 and 9. The volume 
of high-density cold star forming disks forming in DM halos 
was refined to / max = 9. The physical size of mesh cells in the 
simulations was Ax; = 26.16 [10/(1 +z)]2 9 ~ l pc, where / is the 
cell's level of refinement. Each refinement level was integrated 
with its own time step Af/ = Afo2"' w2x 10 4 2 9 "' yrs, where 
Afo < 10 7 yrs is the global time step on the zeroth level, set 
using the CFL condition. 

The two analyzed simulations are the same in all respects, 
except that one included all of the feedback processes described 
above (energy and metal injection, mass loss, etc.), while the 
other did not. The simulation with no feedback is designated 
as NF throughout this paper. Note that the cooling rates in the 
run with feedback accounted for the local metallicity of the gas, 
while in the run with no feedback the significantly lower zero- 
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FIG. 1 . — Star formation rate surface density vs. gas surface density in the 
simulated galaxies at different epochs. The open squares show z = 5 epoch 
in the simulation in which no feedback (NF) processes were included. The 
rates and gas density are averaged over 3.35 [(1 +z)/ 10] kpc. The solid line is 
relation SFR oc S' 4 while the dashed line is SFR oc Eg. 

metallicity cooling rates were used. 

4. RESULTS 

Figure 1 shows the star formation rate per area vs. gas sur- 
face density at different epochs in the simulations. Individual 
points correspond to the level / = 2 cells of size Ax2 = 3.35 [(1 + 
z)/10] kpc. Note that the cells represent not a single object but 
several dozen galaxies with masses > 10 8 /i -1 M Q distributed in 
a high resolution ~ 1 Mpc (comoving) lagrangian region. The 
gas surface density was computed simply as S gas = p g Ax2 and 
the SFR in each cell was estimated as m*t~ x Ax^ 2 , where m* is 
the mass of stars younger than f* = 3 x 10 7 yrs. The averaging 
scale and f* are close to the values typical for observational es- 
timates (see Kennicutt 1998a). The resulting SFR per area and 
E gas are not sensitive to variation in averaging scale and f* over 
a large range of values. 

The figure shows that the star formation rate is strongly cor- 
related with gas surface density: SFR oc X" with n s=s 1 .4 (solid 
line) for £ gas > 10 M Q pc" 2 . At smaller densities the rate rapidly 
decreases and most / = 2 cells with £ gas < 5 M pc~ 2 do not 
contain young stars (i.e., SFR = 0). The simulations thus repro- 
duce the slope of the empirical Schmidt law and the observed 
critical threshold for star formation remarkably well. Figure 1 
also shows that the relation is universal: it does not depend 
sensitively on redshift or input physics. The correlation is not 
sensitive to the variations of averaging scale and f* by at least a 
factor of four. This is a non-trivial result because the correlation 
is observed on kiloparsec scales, while star formation occurs in 
individual 20-50 pc cells with the rate simply proportional to 
the gas density: p\ oc p g . The SFR oc £ g correlation is indeed 
recovered when averaging is done on the scale of individual star 
forming regions (~ 100 pc). The observed scaling SFR oc Eg a | 
on kpc scales is therefore reproduced naturally when a constant 
efficiency of star formation is assumed on the scale of molecu- 



FlG. 2. — Probability distribution function of gas density in the simulated 
gas disks at different redshifts. The PDF, the fraction of total volume in a 
given density interval (Alog| p g = 0.1), was computed in the regions refined 
to / > 8 refinement level. The PDFs are normalized to unity. The histograms 
show outputs at different epochs; the blue histogram shows results for z - 5 
output of the run without feedback (NF). The dashed line is oc pz [ , while the 
dotted line shows a log-normal PDF 

lar complexes. 

The scatter about the mean of the correlation is rather small 
(~ 0.2-0.3 dex) and is well below the observed scatter. The 
observed individual rates and densities, however, have errors 
> 0.2-0.3 dex so larger scatter can be expected. The scatter at 
any given epoch is even smaller, which is remarkable given that 
the galaxies at z > 4 undergo very frequent and violent merg- 
ers. The normalization of the correlation is somewhat lower 
than observed. This can be fixed by lowering r*. The gas con- 
sumption time is a free parameter of the model and so is the 
normalization of the SFR - S g correlation. 

5. DISCUSSION AND CONCLUSIONS 

What mechanism gives rise to the Schmidt law in simulated 
galaxies and what explains its universality? Some understand- 
ing can be gained by considering the density distribution of gas 
in simulated galaxy disks. Figure 2 shows the probability dis- 
tribution function (PDF) of gas density in the disks at different 
epochs. Although the PDF evolves with time, its general shape 
at every epoch can be characterized by the relatively flat region 
at p g < 1 — 10 MqPC -3 and a power law distribution at high 
densities. The high-density tail is also well approximated by 
the log-normal distribution (dotted line in Fig. 2). The evolu- 
tion is driven primarily by shocks associated with supersonic 
turbulence induced by gravitational and shearing instabilities 
and, occasionally, by tidal interactions during mergers in the 
gas disks. For instance, the most massive disk in the simulation 
at z = 4 has well-developed large-scale spiral arms and exhibits 
signs of two minor mergers in progress (Rravtsov 2003). 

Stellar feedback appears to have only a minor effect on the 
overall shape of the PDF. Comparison of PDFs for the feed- 
back and no-feedback (NF) runs at z = 5 shows that distribution 
is remarkably insensitive to the inclusion of several feedback 



4 



processes. The only difference is at the lowest gas densities, 
where the simulation with feedback has considerably more low- 
density hot gas due to the energy injection by SNe. Note that 
in terms of morphology of gas distribution the effect of energy 
injection is significant as the lowest density regions occupy a 
large fraction of the volume. Globally, the effect of energy 
feedback on the thermal state of gas in the halos and surround- 
ing intergalactic medium is also significant (Tassis et al. 2002). 
This makes the insensitivity of the high-density PDF to feed- 
back even more remarkable. The distributions are also not sen- 
sitive to the differences in metallicity and cooling rates between 
feedback and no feedback runs. 

The characteristic shape of the PDF in figure 2 has been 
found in several numerical studies of turbulent ISM. In partic- 
ular, it is very similar to the density PDF in the 2D disk sim- 
ulations of Wada & Norman (2001). These authors also found 
that density PDF is not sensitive to the stellar energy feedback. 
The log-normal high-density tail of the PDF is thought to be a 
generic feature of the turbulent (i.e., randomly forced) super- 
sonic flows (e.g., Vazquez-Semadeni 1994; Padoan et al. 1997; 
Vazquez-Semadeni et al. 2000). The effective pressure in such 
flows is dominated by turbulent (kinematic) rather than ther- 
mal pressure. This may explain the relatively low sensitivity of 
the density distribution to the details of cooling. The cooling 
and heating, however, are important in determining the effec- 
tive equation of state of the gas and thus should definitely also 
play a role. The gas will become increasingly more compress- 
ible as the ratio of the local cooling time to the characteristic 
crossing time decreases. Higher compressibility would then al- 
low the gas to reach higher densities in converging flows, pre- 
sumably creating the high-density tail of the distribution and 
star forming regions. The high-density tail is thus likely to be 
the result of the driven turbulent cascade in a nearly isothermal 
gas. Indeed, for the heating/cooling rates adopted in the simula- 
tions temperatures of the gas at densities p g > 5 M Q pc~ 3 are in 
the range ~ 100- 1000 K. The transition to isothermality, how- 
ever, does not explain the shape of the PDF. The same shape 
was found in a test simulation in which gas was artificially kept 
isothermal at T = 300 K. 

The Schmidt law on kiloparsec scales implies that the ef- 
ficiency of the turbulent build-up of the high-density tail de- 
pends on the average gas density at these scales. Clearly, if 
the PDF shape were independent of average density we would 
have SFR oc £ g . To satisfy the observed non-linear scaling, 
the fraction of mass in star forming regions should scale as 
msp/rrig oc Eg -1 , where n is the slope of the Schmidt law, m g 
is the total mass of gas in a volume element used in averaging. 

Figure 3 shows the density PDFs in 3.35 kpc (/ = 3) cells with 
different densities. Although the low-density part of the distri- 
bution fluctuates, the largest changes are in the high-density tail 
which steepens with decreasing gas density. For £„ = 18 M Q pc~" 
the maximum density is close to the star formation threshold 
Psf = 1-64 MqPC -3 adopted in simulations. The inset in the 
figure 3 shows that the fraction of mass at densities p > p$p 
as a function of average surface density scales as oc £ g 4 , as ex- 
pected from the Schmidt law. Interestingly, this can also explain 
the results of Wong & Blitz (2002) who found that SFR oc S^" 1 
with n mo \ w 1 while SFR oc Shi+h, With n~ 1 .4, where Eh 2 is 
the density of molecular hydrogen. Most of the gas at pressures 
> 10 4 cm" 3 K in their observations is in molecular form. This 
corresponds to densities of > 1-5 M pc~ 3 . The molecular gas 
thus simply traces the high-density tail of the PDF. Figure 3 
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FIG. 3. — Density PDFs in 3.35 kpc (I = 3) mesh cells of different surface 
densities at z = 4. The PDFs are normalized to unity. The dashed line is oc p~' 
dependence. The inset plot shows the fraction of mass at densities p > psp 
as a function of average surface density of cells. The dashed line here shows 
scaling oc Ejj , expected from the Schmidt law (see text for details). 

shows that constant efficiency of the star formation in the dens- 
est regions does not contradict the Schmidt law for the total gas 
density on large scales. 

The scatter in tosf for E g < 10-20 M Q pc~ 2 increases as the 
maximum density of the PDF fluctuates around psf reflecting 
the stochastic nature of turbulent flows. The build-up of the 
high-density tail can be viewed as a random walk of gas ele- 
ments in density with the probability proportional to the density 
PDF (Elmegreen 2002). The bottleneck in the rate of gas cross- 
ing psF is at low densities at largest scales, which therefore con- 
trol the mass of gas available for star formation. Indeed, sim- 
ulations of Wada & Norman (2001) show that the density PDF 
reaches its "equilibrium" shape on the dynamical time corre- 
sponding to the average density of the region. 

This scenario can naturally explain the tight connection be- 
tween the average gas density on kiloparsec scales and star for- 
mation on small scales. The universality of the connection can 
be attributed to the relatively low sensitivity of the density PDF 
to the details of feedback and precise cooling and heating rates. 
The efficiency with which driven supersonic turbulence builds 
up the high-density tail of the PDF depends mainly on the aver- 
age gas density of the region. Results of this study indicate that 
the Schmidt law is simply a manifestation of this dependency. 
It is intriguing that the turbulence may also be responsible for 
the redistribution of angular momentum and exponential light 
profiles of galactic disks (e.g., Lin & Pringle 1987; Silk 2001; 
Slyz et al. 2002). 

The modeling of the density dependence and evolution of the 
PDF is well within the capabilities of the current state-of-the- 
art numerical simulations of the ISM (e.g., Vazquez-Semadeni 
et al. 2000; Wada & Norman 2001; Ostriker et al. 2001). This 
opens the possibility to construct a sound detailed theoretical 
model of the empirical global star formation law in the very 
near future. 
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